%% QR方法---------------------------- 
% clc; clear; close all;
% 
% alpha = 2.4;
% A = [-alpha, 0, 0;
%       0,     0, 0;
%       0,     0, 0];
% B = [1, -4, -3.5;
%      0,  1,  2;
%     -1, -4,  1.5];
% 
% k_list = linspace(0.0001, 0.2, 50);
% dt = 0.01;
% t_transient = 450;
% t_total = 500;
% 
% init_state = [1e-6; 0; 0; 0; 0; 0; 0; 0; 0];
% nLE = 4;
% dim = length(init_state);
% 
% LEs_mat = zeros(length(k_list), nLE);
% 
% parfor idx = 1:length(k_list)
%     k = k_list(idx);
%     LEs_mat(idx, :) = qrLEs(k, A, B, init_state, dt, t_transient, t_total, nLE);
%     fprintf('k=%.3f, LEs= %s\n', k, mat2str(LEs_mat(idx,:),4));
% end
% 
% figure; hold on;
% colors = ['r','b','g','m'];
% for i = 1:nLE
%     plot(k_list, LEs_mat(:,i), 'Color', colors(i), 'LineWidth', 1.5);
% end
% xlabel('k');
% ylabel('LEs');
% legend({'LE_1','LE_2','LE_3','LE_4'});
% title('QR方法');
% grid on; hold off;
% 
% function LEs = qrLEs(k, A, B, x0, dt, t_trans, t_end, nLE)
%     dim = length(x0);
%     h = dt;
%     steps = round(t_end/h);
%     trans_steps = round(t_trans/h);
% 
%     x = x0;
%     Q = eye(dim, nLE); % 初始扰动向量正交矩阵
%     sum_ln = zeros(nLE,1);
%     count = 0;
% 
%     for i = 1:steps
%         % 主轨迹积分
%         x = RK4(@(t,x) mCNN(t,x,A,B,k), x, h);
% 
%         % 计算雅可比矩阵
%         J = Jacobian(@(x)mCNN(0,x,A,B,k), x, 1e-8);
% 
%         % 扰动向量演化
%         V = Q + h * (J * Q);
% 
%         % QR分解正交化
%         [Q,R] = qr(V,0);
% 
%         if i > trans_steps
%             sum_ln = sum_ln + log(abs(diag(R)));
%             count = count + 1;
%         end
%     end
% 
%     LEs = sum_ln / (count * h);
% end


%% Wolf------------------
% clc; clear; close all;
%
% alpha = 2.4;
% A = [-alpha, 0, 0;
%       0,     0, 0;
%       0,     0, 0];
% B = [1, -4, -3.5;
%      0,  1,  2;
%     -1, -4,  1.5];
% 
% k_list = linspace(0.0001, 0.2, 50);
% dt = 0.01;
% t_transient = 450;
% t_total = 500;
% 
% init_state = [1e-6; 0; 0; 0; 0; 0; 0; 0; 0];
% nLE = 4;
% dim = length(init_state);
% 
% LEs_mat = zeros(length(k_list), nLE);
% 
% parfor idx = 1:length(k_list)
%     k = k_list(idx);
%     LEs_mat(idx, :) = wolfLEs_pure(k, A, B, init_state, dt, t_transient, t_total, nLE);
%     fprintf('k=%.3f, LEs= %s\n', k, mat2str(LEs_mat(idx,:),4));
% end
% 
% figure; hold on;
% colors = ['r','b','g','m'];
% for i = 1:nLE
%     plot(k_list, LEs_mat(:,i), 'Color', colors(i), 'LineWidth', 1.5);
% end
% xlabel('k');
% ylabel('LEs');
% legend({'LE_1','LE_2','LE_3','LE_4'});
% title('Wolf方法');
% grid on; hold off;
% 
% function LEs = wolfLEs_pure(k, A, B, x0, dt, t_trans, t_end, nLE)
%     dim = length(x0);
%     h = dt;
%     steps = round(t_end/h);
%     trans_steps = round(t_trans/h);
% 
%     x = x0;
%     Q = eye(dim, nLE); % 多个扰动向量组
%     sum_ln = zeros(nLE,1);
%     count = 0;
% 
%     for i = 1:steps
%         % 主轨迹积分（RK4）
%         x = RK4(@(t,x) mCNN(t,x,A,B,k), x, h);
% 
%         % 计算雅可比矩阵
%         J = Jacobian(@(x)mCNN(0,x,A,B,k), x, 1e-8);
% 
%         % 扰动向量线性演化
%         V = Q + h * (J * Q);
% 
%         % 逐向量归一化并累积增长率（无正交）
%         for j = 1:nLE
%             normV = norm(V(:,j));
%             if normV == 0
%                 normV = 1e-16;
%             end
%             V(:,j) = V(:,j) / normV;
% 
%             if i > trans_steps
%                 sum_ln(j) = sum_ln(j) + log(normV);
%             end
%         end
% 
%         Q = V;
% 
%         if i > trans_steps
%             count = count + 1;
%         end
%     end
% 
%     LEs = sum_ln / (count * h);
% end


%% Gram-Schmidt方法------------------
clc; clear; close all;

% 系统参数
alpha = 2.4;
A = [-alpha, 0, 0;
      0,     0, 0;
      0,     0, 0];
B = [1, -4, -3.5;
     0,  1,  2;
    -1, -4,  1.5];

k_list = linspace(0.0001, 0.2, 50);
dt = 0.01;
t_transient = 200;
t_total = 300;

init_state = [1e-6; 0; 0; 0; 0; 0; 0; 0; 0];
nLE = 4;
dim = length(init_state);

LEs_mat = zeros(length(k_list), nLE);

parfor idx = 1:length(k_list)
    k = k_list(idx);
    LEs_mat(idx,:) = lyapunov_gram_schmidt(k, A, B, init_state, dt, t_transient, t_total, nLE);
    fprintf('k=%.4f, LEs = %s\n', k, mat2str(LEs_mat(idx,:),4));
end

% 绘制
figure; hold on;
colors = ['r','b','g','m'];
for i = 1:nLE
    plot(k_list, LEs_mat(:,i), 'Color', colors(i), 'LineWidth', 1.5);
end
xlabel('k');
ylabel('LEs');
legend({'LE_1','LE_2','LE_3','LE_4'}, 'Location','best');
title('Gram-Schmidt方法');
grid on; hold off;

%% Gram-Schmidt方法计算李指数函数
function LEs = lyapunov_gram_schmidt(k, A, B, x0, dt, t_trans, t_end, nLE)
    dim = length(x0);
    h = dt;
    steps = round(t_end/h);
    trans_steps = round(t_trans/h);

    x = x0;
    Q = eye(dim, nLE);
    sum_ln = zeros(nLE,1);
    count = 0;

    for i = 1:steps
        % RK4积分主轨迹
        x = RK4(@(t,x) mCNN(t,x,A,B,k), x, h);

        % 数值计算雅可比矩阵
        J = Jacobian(@(x)mCNN(0,x,A,B,k), x, 1e-8);

        % 扰动向量演化
        V = Q + h * (J * Q);

        % Gram-Schmidt正交化
        [Q, R] = gram_schmidt(V);

        % 统计对数增长
        if i > trans_steps
            sum_ln = sum_ln + log(abs(diag(R)));
            count = count + 1;
        end
    end

    LEs = sum_ln / (count * h);
end

%% Gram-Schmidt正交化函数（返回正交矩阵Q和上三角矩阵R）
function [Q, R] = gram_schmidt(V)
    [n, m] = size(V);
    Q = zeros(n,m);
    R = zeros(m,m);
    for j = 1:m
        v = V(:,j);
        for i = 1:j-1
            R(i,j) = Q(:,i)' * v;
            v = v - R(i,j)*Q(:,i);
        end
        R(j,j) = norm(v);
        Q(:,j) = v / R(j,j);
    end
end
